library(foreign)
library(readr)
library(survival)
library(haven)
library(MASS)
library(modeLLtest)
library(coxed)
library(ggplot2)
library(coxrobust)
library(stargazer)
library(xtable)
library(lme4)
library(plm)
library(gplots)
library(arm)
library(knitr)
library(nnet)
library(gmodels)
library(expss)
library(dplyr) 
library(tidyr)
library(xtable)
library(reporttools)
library(survey)
library(marginaleffects)
library(devtools)
library(ggeffects)
setwd("~/Desktop")

##################################
##########Loading Data############ 
##################################

data<-read_csv("Dataset_Early_Bagdanov.csv")


####################################################
########Percentages of Individuals who Report#######  
#############Engaging (Table 3)#####################
####################################################

#Transport- Buses

count_buses<-length(which(data$isr.buses==1))
percent_buses<-(count_buses/1255)*100
percent_buses

#Transport- Lightrail 

count_lightrail<-length(which(data$isr.lightrail==1))
percent_lightrail<-(count_lightrail/1255)*100
percent_lightrail

#Police

count_police<-length(which(data$police==1))
percent_police<-(count_police/1255)*100
percent_police

#Police Station

count_police_station<-length(which(data$police_station==1))
percent_police_station<-(count_police_station/1255)*100
percent_police_station

#Dispute Res

count_dispute<-length(which(data$dispute==1))
percent_dispute<-(count_dispute/1255)*100
percent_dispute

#Justice_crime

count_justice_crime<-length(which(data$justice_crime==1))
percent_justice_crime<-(count_justice_crime/1255)*100
percent_justice_crime

#Justice_victim

count_justice_victim<-length(which(data$justice_victim==1))
percent_justice_victim<-(count_justice_victim/1255)*100
percent_justice_victim

#School_Self 

count_school_self<-length(which(data$bagrut==1))
percent_school_self<-(count_school_self/1255)*100
percent_school_self

#School_Children

count_edu_child<-length(which(data$edu_child==1))
percent_edu_child<-(count_edu_child/813)*100
percent_edu_child

#Higher ed 

count_hu<-length(which(data$HU==1))
percent_hu<-(count_hu/1255)*100
percent_hu

#Healthcare (Primary Care Clinic) 

count_clinic<-length(which(data$clinic==1))
percent_clinic<-(count_clinic/1255)*100
percent_clinic

#Health (Prenatal Care) 

count_prenatal_hospital<-length(which(data$hospital==1))
percent_prenatal_hospital<-(count_prenatal_hospital/1255)*100
percent_prenatal_hospital

#Economic Assistance 

count_econ_assist<-length(which(data$econ_assist==1))
percent_econ_assist<-(count_econ_assist/1255)*100
percent_econ_assist


#Housing 

count_housing<-length(which(data$housing==1))
percent_housing<-(count_housing/1255)*100
percent_housing

#Water 

count_water<-length(which(data$water==1))
percent_water<-(count_water/1255)*100
percent_water

#Sanitation 

count_sanitation<-length(which(data$sanitation==1))
percent_sanitation<-(count_sanitation/1255)*100
percent_sanitation

#Infrastructure

count_infrastructure<-length(which(data$infrastructure==1))
percent_infrastructure<-(count_infrastructure/1255)*100
percent_infrastructure

#Electricity 

count_electricity<-length(which(data$electricity==1))
percent_electricity<-(count_electricity/1255)*100
percent_electricity

#Community Center

count_cc<-length(which(data$comm_center==1))
percent_cc<-(count_cc/1255)*100
percent_cc

#Parks

count_parks<-length(which(data$parks==1))
percent_parks<-(count_parks/1255)*100
percent_parks

#Voting

count_voting<-length(which(data$voting==1))
percent_voting<-(count_voting/1255)*100
percent_voting

####################################################
########Logit Regressions (Tables 4 and 5)##########
####################################################

#Police

##With al_fadi 
police.glm<-glm(police_station ~ police_norm_bi + state_respon + sex + as.numeric(age) + region + edu_level + al_fadi+ married + relig_attend + income + refugee + prison + party, data=data, binomial(link = "logit")) 
summary(police.glm)

##Without al_fadi for appendix
m1<-glm(police_station ~ police_norm_bi + state_respon + sex + as.numeric(age) + region + edu_level + married + relig_attend + income + refugee + prison + party, data=data, binomial(link = "logit")) 
summary(m1)

#Transport 

##With al_fadi
transport.glm<-glm(transport~ transport_norm_bi + state_respon +sex + as.numeric(age) + region + edu_level + al_fadi+ married + relig_attend + income + refugee + prison + party, data=data, binomial(link = "logit")) 
summary(transport.glm)

##Without al_fadi for appendix
m2<-glm(transport~ transport_norm_bi + state_respon +sex + as.numeric(age) + region + edu_level +  married + relig_attend + income + refugee + prison + party, data=data, binomial(link = "logit")) 
summary(m2)

#Education
##With al_fadi
edu.glm<-glm(edu_isr~ edu_norm_bi + state_respon + sex + as.numeric(age) + region + edu_level + al_fadi + married + relig_attend + income + refugee + prison + party, data=data, binomial(link = "logit")) 
summary(edu.glm)

##Without al_fadi for appendix
m3<-glm(edu_isr~ edu_norm_bi + state_respon + sex + as.numeric(age) + region + edu_level + married + relig_attend + income + refugee + prison + party, data=data, binomial(link = "logit")) 
summary(m3)

#Healthcare

##With al_fadi
medical.glm<-glm(clinic~ medical_norm_bi + state_respon + sex + as.numeric(age) + region + edu_level + al_fadi+ married + relig_attend + income + refugee + prison + party, data=data, binomial(link = "logit")) 
summary(medical.glm)

##Without al_fadi for appendix
m4<-glm(clinic~ medical_norm_bi + state_respon + sex + as.numeric(age) + region + edu_level +  married + relig_attend + income + refugee + prison + party, data=data, binomial(link = "logit")) 
summary(m4)

#Community Center

##With al_fadi
cc.glm<-glm(comm_center~ cc_norm_bi + state_respon+ sex + as.numeric(age) + region + edu_level + al_fadi+ married + relig_attend + income + refugee + prison + party, data=data, binomial(link = "logit")) 
summary(cc.glm)

##Without al_fadi for appendix
m5<-glm(comm_center~ cc_norm_bi + state_respon+ sex + as.numeric(age) + region + edu_level + married + relig_attend + income + refugee + prison + party, data=data, binomial(link = "logit")) 
summary(m5)

#Sanitation

##With al_fadi
sanitation.glm<-glm(sanitation~ sanitation_norm_bi + state_respon+ sex + as.numeric(age) + region + edu_level + al_fadi+ married + relig_attend + income + refugee + prison + party, data=data, binomial(link = "logit")) 
summary(sanitation.glm)

##Without al_fadi for appendix
m6<-glm(sanitation~ sanitation_norm_bi + state_respon+ sex + as.numeric(age) + region + edu_level + married + relig_attend + income + refugee + prison + party, data=data, binomial(link = "logit")) 
summary(m6)

#Voting

##With al_fadi
voting.glm<-glm(voting ~ voting_norm_bi +state_respon+ sex + as.numeric(age) + region + edu_level + al_fadi+ married + relig_attend + income + refugee + prison + party, data=data, binomial(link = "logit")) 
summary(voting.glm)

##Without al_fadi for appendix
m7<-glm(voting ~ voting_norm_bi +state_respon+ sex + as.numeric(age) + region + edu_level +  married + relig_attend + income + refugee + prison + party, data=data, binomial(link = "logit")) 
summary(m7)

#Welfare/Economic Assistance

##With al_fadi
econ.glm<-glm(econ_assist~ econ_assist_norm_bi+ state_respon+sex + as.numeric(age) + region + edu_level + al_fadi+ married + relig_attend + income + refugee + prison + party, data=data, binomial(link = "logit")) 
summary(econ.glm)

##Without al_fadi for appendix
m8<-glm(econ_assist~ econ_assist_norm_bi+ state_respon+sex + as.numeric(age) + region + edu_level + married + relig_attend + income + refugee + prison + party, data=data, binomial(link = "logit")) 
summary(m8)

#Parks 

##With al_fadi
park.glm<-glm(parks~ park_norm_bi+state_respon+sex + as.numeric(age)  + region + edu_level + al_fadi+ married + relig_attend + income + refugee + prison + party, data=data, binomial(link = "logit")) 
summary(park.glm)

##Without al_fadi for appendix
m9<-glm(parks~ park_norm_bi+state_respon+sex + as.numeric(age)  + region + edu_level + married + relig_attend + income + refugee + prison + party, data=data, binomial(link = "logit")) 
summary(m9)

#Dispute Res

##With al_fadi
dispute.glm<-glm(dispute~dispute_norm_bi + state_respon+ sex + as.numeric(age) + region + edu_level + al_fadi+ married + relig_attend + income + refugee + prison + party, data=data, binomial(link = "logit")) 
summary(dispute.glm)

##Without al_fadi for appendix
m10<-glm(dispute~dispute_norm_bi + state_respon+ sex + as.numeric(age) + region + edu_level + married + relig_attend + income + refugee + prison + party, data=data, binomial(link = "logit")) 
summary(m10)

#Generating Latex Tables (Table 4 and Table 5)

stargazer(police.glm,edu.glm,voting.glm, sanitation.glm,econ.glm, star.cutoffs = c(0.05, 0.01, 0.001))
stargazer(park.glm,dispute.glm,cc.glm,medical.glm,transport.glm,star.cutoffs = c(0.05, 0.01, 0.001))

#Generating Latex table for Supplementary Materials Section F (Table 8)
stargazer(m1,m2,m3,m4,m5,m6,m7,m8,m9,m10, star.cutoffs = c(0.05, 0.01, 0.001))

#################################################
#############Predicted Probabilities#############
#################(Figure 3)######################
#################################################

#Police

##Make IV into a factor 
data$police_fac<-as.factor(data$police_norm_bi)

##Rerun GLM with factor variable
police.glm.fac<-glm(police_station~ police_fac + state_respon + sex  + region + edu_level + al_fadi + married + as.numeric(age) + relig_attend + income + refugee + prison + party, data=data, binomial(link = "logit")) 
summary(police.glm.fac)

##Generate Predicted Probabilities 
ggpredict(police.glm.fac,terms = "police_fac")

##Generate Plot 
police.plot<-ggpredict(police.glm.fac,terms = "police_fac") %>%
  ggplot(mapping=aes(x=x,y=predicted))+
  geom_point()+
  labs(title="Visiting Police Station",x="Normalization Perception",y="Predicted Probability")+
  theme(plot.title = element_text(hjust = 0.5, size=12))+
  scale_y_continuous(limits = c(0.10, 1),breaks = seq(0.10, 1, by = .05)) +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=.5,alpha=.3)+geom_segment(x=1,y=.35,xend=2, yend=0.29)
police.plot

#Education

##Make IV into a factor 
data$edu_fac<-as.factor(data$edu_norm_bi)

##Rerun GLM with factor variable
edu.glm.fac<-glm(edu_isr~ edu_fac +state_respon + sex  + region + edu_level +as.numeric(age)+al_fadi  +married + relig_attend + income + refugee + prison + party, data=data, binomial(link = "logit")) 
summary(edu.glm.fac)

##Generate Predicted Probabilities 
ggpredict(edu.glm.fac,terms = "edu_fac")

##Generate Plot  
edu.plot<-ggpredict(edu.glm.fac,terms = "edu_fac") %>%
  ggplot(mapping=aes(x=x,y=predicted))+
  geom_point()+ geom_line()+
  labs(title="Attending Jerusalem Municipality School",x="Normalization Perception",y="Predicted Probability")+
  geom_smooth(se=FALSE) +
  theme(plot.title = element_text(hjust = 0.5, size=12))+
  scale_y_continuous(limits = c(0.10, 1),breaks = seq(0.10, 1, by = .05)) +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=.5,alpha=.3)+geom_segment(x=1,y=.63,xend=2, yend=0.51)
edu.plot

#Sanitation

##Make IV into a factor 
data$sanitation_fac<-as.factor(data$sanitation_norm_bi)

##Rerun GLM with factor variable
sanitation.glm.fac<-glm(sanitation~ sanitation_fac + state_respon+ sex  + region + edu_level + as.numeric(age) + al_fadi+ married + relig_attend + income + refugee + prison + party, data=data, binomial(link = "logit")) 
summary(sanitation.glm.fac)

##Generate Predicted Probabilities 
ggpredict(sanitation.glm.fac,terms = "sanitation_fac")


##Generate Plot
sanitation.plot<-ggpredict(sanitation.glm.fac,terms = "sanitation_fac") %>%
  ggplot(mapping=aes(x=x,y=predicted))+
  geom_point()+
  labs(title="Contacting Sanitation Department",x="Normalization Perception",y="Predicted Probability")+
  geom_smooth(se=FALSE) + 
  theme(plot.title = element_text(hjust = 0.5, size=12))+
  scale_y_continuous(limits = c(0.10, 1),breaks = seq(0.10, 1, by = .05)) +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=.5,alpha=.3)+geom_segment(x=1,y=.96,xend=2, yend=0.88)
sanitation.plot  

#Voting

##Make IV into a factor 
data$voting_fac<-as.factor(data$voting_norm_bi)

##Rerun GLM with factor variable
voting.glm.fac<-glm(voting ~ voting_fac + state_respon+ sex  + region + edu_level +al_fadi + married + relig_attend + as.numeric(age)+ income + refugee + prison + party, data=data, binomial(link = "logit")) 
summary(voting.glm.fac)

##Generate Predicted Probabilities 
ggpredict(voting.glm.fac,terms = "voting_fac")

##Generate Plot
voting.plot<-ggpredict(voting.glm.fac,terms = "voting_fac") %>%
  ggplot(mapping=aes(x=x,y=predicted))+
  geom_point()+
  labs(title="Voting",x="Normalization Perception",y="Predicted Probability")+
  geom_smooth(se=FALSE) + 
  theme(plot.title = element_text(hjust = 0.5, size=12))+
  scale_y_continuous(limits = c(0.10, 1),breaks = seq(0.10, 1, by = .05)) +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=.5,alpha=.3)+geom_segment(x=1,y=.37,xend=2, yend=0.23)
voting.plot 

#Econ Assistance

##Make  IV into a factor 
data$econ_fac<-as.factor(data$econ_assist_norm_bi)

##Rerun GLM with factor variable
econ.glm.fac<-glm(econ_assist~ econ_fac + state_respon+sex  + region + edu_level + al_fadi+ married + relig_attend + income + as.numeric(age) + refugee + prison + party, data=data, binomial(link = "logit")) 
summary(econ.glm.fac)

##Generate Predicted Probabilities 
ggpredict(econ.glm.fac,terms = "econ_fac")

##Generate Plot 
econ.plot<-ggpredict(econ.glm.fac,terms = "econ_fac") %>%
  ggplot(mapping=aes(x=x,y=predicted))+
  geom_point()+
  labs(title="Seeking \nSocial Welfare/Economic Assistance",x="Normalization Perception",y="Predicted Probability")+
  geom_smooth(se=FALSE) + 
  theme(plot.title = element_text(hjust = 0.5, size=12))+
  scale_y_continuous(limits = c(0.10, 1),breaks = seq(0.10, 1, by = .05))+
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=.5,alpha=.3)+geom_segment(x=1,y=.62,xend=2, yend=0.24)
econ.plot 

##Arrange in grid to generate Figure 3
library(gridExtra)
grid.arrange(econ.plot, voting.plot, police.plot,sanitation.plot, edu.plot, ncol=3, nrow=2)

#############################################
################Experiment###################
################(Figure 4)###################
#############################################

#Loading Experiment Data with Question Iterations stacked in single column (7530 observations)

merge_final<-read_csv("exp_replicate.csv")
View(merge_final)

#Counting likeability, responsibility, respectibility, and neighborliness scores and creating composite variable

##Likability 

###Police
merge_final$plike<-rep(NA,length(7530))
merge_final$plike[merge_final$Numeric_1==3]<-merge_final$Q29a_1 #[merge_final$Numeric_1==3]
merge_final$plike[merge_final$Numeric_2==3]<-merge_final$Q29b_1 #[merge_final$Numeric_2==3]
merge_final$plike[merge_final$Numeric_3==3]<-merge_final$Q29c_1 #[merge_final$Numeric_3==3]
merge_final$plike[merge_final$Numeric_4==3]<-merge_final$Q29d_1 #[merge_final$Numeric_4==3]
merge_final$plike[merge_final$Numeric_5==3]<-merge_final$Q29e_1 #[merge_final$Numeric_5==3]
merge_final$plike[merge_final$Numeric_6==3]<-merge_final$Q29f_1 #[merge_final$Numeric_6==3]
table(merge_final$plike)

###Health 
merge_final$hlike<-rep(NA,length(7530))
merge_final$hlike[merge_final$Numeric_1==1]<-merge_final$Q29a_1[merge_final$Numeric_1==1]
merge_final$hlike[merge_final$Numeric_2==1]<-merge_final$Q29b_1[merge_final$Numeric_2==1]
merge_final$hlike[merge_final$Numeric_3==1]<-merge_final$Q29c_1[merge_final$Numeric_3==1]
merge_final$hlike[merge_final$Numeric_4==1]<-merge_final$Q29d_1[merge_final$Numeric_4==1]
merge_final$hlike[merge_final$Numeric_5==1]<-merge_final$Q29e_1[merge_final$Numeric_5==1]
merge_final$hlike[merge_final$Numeric_6==1]<-merge_final$Q29f_1[merge_final$Numeric_6==1]
table(merge_final$hlike)

###Vote 
merge_final$vlike<-rep(NA,length(7530))
merge_final$vlike[merge_final$Numeric_1==4]<-merge_final$Q29a_1[merge_final$Numeric_1==4]
merge_final$vlike[merge_final$Numeric_2==4]<-merge_final$Q29b_1[merge_final$Numeric_2==4]
merge_final$vlike[merge_final$Numeric_3==4]<-merge_final$Q29c_1[merge_final$Numeric_3==4]
merge_final$vlike[merge_final$Numeric_4==4]<-merge_final$Q29d_1[merge_final$Numeric_4==4]
merge_final$vlike[merge_final$Numeric_5==4]<-merge_final$Q29e_1[merge_final$Numeric_5==4]
merge_final$vlike[merge_final$Numeric_6==4]<-merge_final$Q29f_1[merge_final$Numeric_6==4]
table(merge_final$vlike)


###Community center 
merge_final$cclike<-rep(NA,length(7530))
merge_final$cclike[merge_final$Numeric_1==2]<-merge_final$Q29a_1[merge_final$Numeric_1==2]
merge_final$cclike[merge_final$Numeric_2==2]<-merge_final$Q29b_1[merge_final$Numeric_2==2]
merge_final$cclike[merge_final$Numeric_3==2]<-merge_final$Q29c_1[merge_final$Numeric_3==2]
merge_final$cclike[merge_final$Numeric_4==2]<-merge_final$Q29d_1[merge_final$Numeric_4==2]
merge_final$cclike[merge_final$Numeric_5==2]<-merge_final$Q29e_1[merge_final$Numeric_5==2]
merge_final$cclike[merge_final$Numeric_6==2]<-merge_final$Q29f_1[merge_final$Numeric_6==2]
table(merge_final$cclike)

##Responsibility

###Police 

merge_final$prespon<-rep(NA,length(7530))
merge_final$prespon[merge_final$Numeric_1==3]<-merge_final$Q29a_2[merge_final$Numeric_1==3]
merge_final$prespon[merge_final$Numeric_2==3]<-merge_final$Q29b_2[merge_final$Numeric_2==3]
merge_final$prespon[merge_final$Numeric_3==3]<-merge_final$Q29c_2[merge_final$Numeric_3==3]
merge_final$prespon[merge_final$Numeric_4==3]<-merge_final$Q29d_2[merge_final$Numeric_4==3]
merge_final$prespon[merge_final$Numeric_5==3]<-merge_final$Q29e_2[merge_final$Numeric_5==3]
merge_final$prespon[merge_final$Numeric_6==3]<-merge_final$Q29f_2[merge_final$Numeric_6==3]
table(merge_final$prespon)

##Health 
merge_final$hrespon<-rep(NA,length(7530))
merge_final$hrespon[merge_final$Numeric_1==1]<-merge_final$Q29a_2[merge_final$Numeric_1==1]
merge_final$hrespon[merge_final$Numeric_2==1]<-merge_final$Q29b_2[merge_final$Numeric_2==1]
merge_final$hrespon[merge_final$Numeric_3==1]<-merge_final$Q29c_2[merge_final$Numeric_3==1]
merge_final$hrespon[merge_final$Numeric_4==1]<-merge_final$Q29d_2[merge_final$Numeric_4==1]
merge_final$hrespon[merge_final$Numeric_5==1]<-merge_final$Q29e_2[merge_final$Numeric_5==1]
merge_final$hrespon[merge_final$Numeric_6==1]<-merge_final$Q29f_2[merge_final$Numeric_6==1]
table(merge_final$hrespon)


##Vote 

merge_final$vrespon<-rep(NA,length(7530))
merge_final$vrespon[merge_final$Numeric_1==4]<-merge_final$Q29a_2[merge_final$Numeric_1==4]
merge_final$vrespon[merge_final$Numeric_2==4]<-merge_final$Q29b_2[merge_final$Numeric_2==4]
merge_final$vrespon[merge_final$Numeric_3==4]<-merge_final$Q29c_2[merge_final$Numeric_3==4]
merge_final$vrespon[merge_final$Numeric_4==4]<-merge_final$Q29d_2[merge_final$Numeric_4==4]
merge_final$vrespon[merge_final$Numeric_5==4]<-merge_final$Q29e_2[merge_final$Numeric_5==4]
merge_final$vrespon[merge_final$Numeric_6==4]<-merge_final$Q29f_2[merge_final$Numeric_6==4]
table(merge_final$vrespon)

##Community Center
merge_final$ccrespon<-rep(NA,length(7530))
merge_final$ccrespon[merge_final$Numeric_1==2]<-merge_final$Q29a_2[merge_final$Numeric_1==2]
merge_final$ccrespon[merge_final$Numeric_2==2]<-merge_final$Q29b_2[merge_final$Numeric_2==2]
merge_final$ccrespon[merge_final$Numeric_3==2]<-merge_final$Q29c_2[merge_final$Numeric_3==2]
merge_final$ccrespon[merge_final$Numeric_4==2]<-merge_final$Q29d_2[merge_final$Numeric_4==2]
merge_final$ccrespon[merge_final$Numeric_5==2]<-merge_final$Q29e_2[merge_final$Numeric_5==2]
merge_final$ccrespon[merge_final$Numeric_6==2]<-merge_final$Q29f_2[merge_final$Numeric_6==2]
table(merge_final$ccrespon)

##Respectability

###Police
merge_final$prespect<-rep(NA,length(7530))
merge_final$prespect[merge_final$Numeric_1==3]<-merge_final$Q29a_3[merge_final$Numeric_1==3]
merge_final$prespect[merge_final$Numeric_2==3]<-merge_final$Q29b_3[merge_final$Numeric_2==3]
merge_final$prespect[merge_final$Numeric_3==3]<-merge_final$Q29c_3[merge_final$Numeric_3==3]
merge_final$prespect[merge_final$Numeric_4==3]<-merge_final$Q29d_3[merge_final$Numeric_4==3]
merge_final$prespect[merge_final$Numeric_5==3]<-merge_final$Q29e_3[merge_final$Numeric_5==3]
merge_final$prespect[merge_final$Numeric_6==3]<-merge_final$Q29f_3[merge_final$Numeric_6==3]
table(merge_final$prespect)

###Health 

merge_final$hrespect<-rep(NA,length(7530))
merge_final$hrespect[merge_final$Numeric_1==1]<-merge_final$Q29a_3[merge_final$Numeric_1==1]
merge_final$hrespect[merge_final$Numeric_2==1]<-merge_final$Q29b_3[merge_final$Numeric_2==1]
merge_final$hrespect[merge_final$Numeric_3==1]<-merge_final$Q29c_3[merge_final$Numeric_3==1]
merge_final$hrespect[merge_final$Numeric_4==1]<-merge_final$Q29d_3[merge_final$Numeric_4==1]
merge_final$hrespect[merge_final$Numeric_5==1]<-merge_final$Q29e_3[merge_final$Numeric_5==1]
merge_final$hrespect[merge_final$Numeric_6==1]<-merge_final$Q29f_3[merge_final$Numeric_6==1]
table(merge_final$hrespect)

###Vote 

merge_final$vrespect<-rep(NA,length(7530))
merge_final$vrespect[merge_final$Numeric_1==4]<-merge_final$Q29a_3[merge_final$Numeric_1==4]
merge_final$vrespect[merge_final$Numeric_2==4]<-merge_final$Q29b_3[merge_final$Numeric_2==4]
merge_final$vrespect[merge_final$Numeric_3==4]<-merge_final$Q29c_3[merge_final$Numeric_3==4]
merge_final$vrespect[merge_final$Numeric_4==4]<-merge_final$Q29d_3[merge_final$Numeric_4==4]
merge_final$vrespect[merge_final$Numeric_5==4]<-merge_final$Q29e_3[merge_final$Numeric_5==4]
merge_final$vrespect[merge_final$Numeric_6==4]<-merge_final$Q29f_3[merge_final$Numeric_6==4]
table(merge_final$vrespect)

###Community Center
merge_final$ccrespect<-rep(NA,length(7530))
merge_final$ccrespect[merge_final$Numeric_1==2]<-merge_final$Q29a_3[merge_final$Numeric_1==2]
merge_final$ccrespect[merge_final$Numeric_2==2]<-merge_final$Q29b_3[merge_final$Numeric_2==2]
merge_final$ccrespect[merge_final$Numeric_3==2]<-merge_final$Q29c_3[merge_final$Numeric_3==2]
merge_final$ccrespect[merge_final$Numeric_4==2]<-merge_final$Q29d_3[merge_final$Numeric_4==2]
merge_final$ccrespect[merge_final$Numeric_5==2]<-merge_final$Q29e_3[merge_final$Numeric_5==2]
merge_final$ccrespect[merge_final$Numeric_6==2]<-merge_final$Q29f_3[merge_final$Numeric_6==2]
table(merge_final$ccrespect)


##Neighborlinesss

###Police
merge_final$pneigh<-rep(NA,length(7530))
merge_final$pneigh[merge_final$Numeric_1==3]<-merge_final$Q29a_4[merge_final$Numeric_1==3]
merge_final$pneigh[merge_final$Numeric_2==3]<-merge_final$Q29b_4[merge_final$Numeric_2==3]
merge_final$pneigh[merge_final$Numeric_3==3]<-merge_final$Q29c_4[merge_final$Numeric_3==3]
merge_final$pneigh[merge_final$Numeric_4==3]<-merge_final$Q29d_4[merge_final$Numeric_4==3]
merge_final$pneigh[merge_final$Numeric_5==3]<-merge_final$Q29e_4[merge_final$Numeric_5==3]
merge_final$pneigh[merge_final$Numeric_6==3]<-merge_final$Q29f_4[merge_final$Numeric_6==3]
table(merge_final$pneigh)

###Health 

merge_final$hneigh<-rep(NA,length(7530))
merge_final$hneigh[merge_final$Numeric_1==1]<-merge_final$Q29a_4[merge_final$Numeric_1==1]
merge_final$hneigh[merge_final$Numeric_2==1]<-merge_final$Q29b_4[merge_final$Numeric_2==1]
merge_final$hneigh[merge_final$Numeric_3==1]<-merge_final$Q29c_4[merge_final$Numeric_3==1]
merge_final$hneigh[merge_final$Numeric_4==1]<-merge_final$Q29d_4[merge_final$Numeric_4==1]
merge_final$hneigh[merge_final$Numeric_5==1]<-merge_final$Q29e_4[merge_final$Numeric_5==1]
merge_final$hneigh[merge_final$Numeric_6==1]<-merge_final$Q29f_4[merge_final$Numeric_6==1]
table(merge_final$hneigh)

###Vote 

merge_final$vneigh<-rep(NA,length(7530))
merge_final$vneigh[merge_final$Numeric_1==4]<-merge_final$Q29a_4[merge_final$Numeric_1==4]
merge_final$vneigh[merge_final$Numeric_2==4]<-merge_final$Q29b_4[merge_final$Numeric_2==4]
merge_final$vneigh[merge_final$Numeric_3==4]<-merge_final$Q29c_4[merge_final$Numeric_3==4]
merge_final$vneigh[merge_final$Numeric_4==4]<-merge_final$Q29d_4[merge_final$Numeric_4==4]
merge_final$vneigh[merge_final$Numeric_5==4]<-merge_final$Q29e_4[merge_final$Numeric_5==4]
merge_final$vneigh[merge_final$Numeric_6==4]<-merge_final$Q29f_4[merge_final$Numeric_6==4]
table(merge_final$vneigh)

###Community Center
merge_final$ccneigh<-rep(NA,length(7530))
merge_final$ccneigh[merge_final$Numeric_1==2]<-merge_final$Q29a_4[merge_final$Numeric_1==2]
merge_final$ccneigh[merge_final$Numeric_2==2]<-merge_final$Q29b_4[merge_final$Numeric_2==2]
merge_final$ccneigh[merge_final$Numeric_3==2]<-merge_final$Q29c_4[merge_final$Numeric_3==2]
merge_final$ccneigh[merge_final$Numeric_4==2]<-merge_final$Q29d_4[merge_final$Numeric_4==2]
merge_final$ccneigh[merge_final$Numeric_5==2]<-merge_final$Q29e_4[merge_final$Numeric_5==2]
merge_final$ccneigh[merge_final$Numeric_6==2]<-merge_final$Q29f_4[merge_final$Numeric_6==2]
table(merge_final$ccneigh)


##Creating Composite Score (likability, respectability, responsibility, neighborliness)

###Police
merge_final$psum_long<-rep(NA,length(7530))
merge_final$psum_long[merge_final$Numeric_1==3]<-merge_final$Q29a_sum[merge_final$Numeric_1==3]
table(merge_final$psum_long)


###Health 
merge_final$hsum_long<-rep(NA,length(7530))
merge_final$hsum_long[merge_final$Numeric_1==1]<-merge_final$Q29a_sum[merge_final$Numeric_1==1]
table(merge_final$hsum_long)

###Vote 
merge_final$vsum_long<-rep(NA,length(7530))
merge_final$vsum_long[merge_final$Numeric_1==4]<-merge_final$Q29a_sum[merge_final$Numeric_1==4]
table(merge_final$vsum_long)


###Community Center
merge_final$ccsum_long<-rep(NA,length(7530))
merge_final$ccsum_long[merge_final$Numeric_1==2]<-merge_final$Q29a_sum[merge_final$Numeric_1==2]
table(merge_final$ccsum_long)

##Creating dataframe 
temp1<-merge_final[,c("hsum_long")]
temp1$attribute<-"health"
temp1<-dplyr::rename(temp1,composite=hsum_long)
temp2<-merge_final[,c("psum_long")]
temp2$attribute<-"police"
temp2<-dplyr::rename(temp2,composite=psum_long)
temp3<-merge_final[,c("vsum_long")]
temp3$attribute<-"vote"
temp3<-dplyr::rename(temp3,composite=vsum_long)
temp4<-merge_final[,c("ccsum_long")]
temp4$attribute<-"community center"
temp4<-dplyr::rename(temp4,composite=ccsum_long)
temp5<-rbind(temp1,temp2,temp3,temp4)
temp6<-na.omit(temp5)
View(temp6)
temp7<-cbind(temp6,merge_final$Name_1,merge_final$Age_1,merge_final$Occupation_1,merge_final$Neighborhood_1,merge_final$SbjNum)
temp7

##Relevel attribute variable so that health is assigned as the control 

temp7$attribute <- factor(temp7$attribute, levels = c("health", "community center", "police","vote"))
temp7

#Running the ATE model
library(sandwich)
library(lmtest)

amce<-lm(composite~attribute+merge_final$Name_1+merge_final$Age_1+merge_final$Occupation_1+merge_final$Neighborhood_1+merge_final$SbjNum, data=temp7, is.na=TRUE)
summary(amce)

##Cluster standard errors 
cluster.se<-coeftest(amce,vcov=vcovCL, cluster=~merge_final$SbjNum)
summary(cluster.se)

##Generate Regression output table (Supplementary Materials Table 3)
stargazer(cluster.se, star.cutoffs = c(0.05, 0.01, 0.001))

##Create Dataframe 
c2 <- confint(cluster.se) #Confidence intervals from regression
cluster.se2 <- as.data.frame(coef(cluster.se)) #put coefficients output into a data frame
cluster.se2  $ci95_b <- c2[,1]# Put confidence intervals into the data frame 
cluster.se2  $ci95_t <- c2[,2]# Put confidence intervals into the data frame 
cluster.se2 $measure <- "composite" #make a measure indicator
cluster.se2 $cat <-"Aggregate"
cluster.se2 $Attribute <- c("Health", "Community Center","Police","Vote","Name 1","Name 2","Name 3","Name 4","Name 5","Name 6","Name 7","Age", "Doctor","Engineer","School Principal","Beit Hanina","Ras Al Amoud","Sur Baher","Sbj_num") #make an attribute indicator
colnames(cluster.se2 )[1] ="mean"
View(cluster.se2)

##Create Plot (Figure 4)

cluster.se2.omit<-cluster.se2[-(5:19),]
cluster.se.omit.2<-cluster.se2.omit[-(1),]
Fig4 <- ggplot(cluster.se.omit.2, aes(x = factor(measure), y = mean)) + 
  geom_point( ) + 
  geom_errorbar( aes(ymin = ci95_b, ymax= ci95_t), width = 0.1) +
  facet_wrap(~Attribute, ncol=4) +
  theme_bw() + #theme(axis.title.x = theme_text(size = 12, vjust = .25))+ 
  scale_x_discrete(breaks=c("composite","composite","composite"),labels=c()) +
  xlab("GSI") +  ylab("Estimated Effect of GSI Usage on Composite Likability")  +scale_y_continuous(limit = c(-9, 0)) # breaks=c(-0.4, -0.2,-0.1, 0, 0.1,0.2,  0.4, 0.6 )) 
Fig4


#################################################
#################Figure 1 Map####################
#################################################

library(sf)
library(RColorBrewer)
library(sp)
library(ggspatial)

#Download the following files from Dataverse and save them to your Desktop, which is the working directory specified on line 28 of this code

##lamas_census_tracts_2011_jerusalem.shx
##lamas_census_tracts_2011_jerusalem.shp
##lamas_census_tracts_2011_jerusalem.dbf
##lamas_census_tracts_2011_jerusalem.prj

#Load Shape File
shape<-read_sf('lamas_census_tracts_2011_jerusalem.shp')

#Load csv with extra measures
map_measures<-read_csv("maps_data.csv")

#Merge files back together 
merged_shp<-sp::merge(shape,map_measures)

#Generate East vs. West Jerusalem Map
ggplot() +
  layer_spatial(
    merged_shp, 
    aes(fill = as.factor((East_West)))
  ) +
  scale_fill_manual(
    breaks = c(1, 0),
    values = c("gray", "white"),
    labels=c("East","West")
  )+ labs(
    x = NULL, 
    y = NULL, 
    title = "East vs. West Jerusalem (1967-Present)", 
    subtitle = "Statistical Areas", 
    caption = "",
    fill = "Region"
  ) + theme(axis.text.x = element_blank(),
            axis.text.y = element_blank(), 
            panel.background= element_rect(fill=("gray90")),
            plot.background = element_rect(fill=("gray90")),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            axis.ticks = element_blank())



